Background

Using 2018-2020 SHIP life expectancy data and 2016-2020 ACS socioeconomic characteristics data of Maryland Counties, I will be conducting a logistic regression in R predicting the odds of a given county having above state average life expectancy, 78.6 years, based on median household income, education level, racial composition, and sex composition.

Mapping the Data

library(tidyverse)  
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)         
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(mapview)    
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)
library(leaflet) 
library(dplyr)
library(leafpop)
library(readr)
md <- read_sf("/Users/mkheyfets/Documents/life_exp/Maryland_Physical_Boundaries_-_County_Boundaries_(Generalized)")
life_exp <- read_csv("scatterplot/SHIP_Life_Expectancy_2010-2020_20231228.csv")
Rows: 24 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Jurisdiction
dbl (1): Life_Expectancy

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
life_exp_j <- life_exp %>% 
  rename("county" = "Jurisdiction") 
life_exp_md <- full_join(md, life_exp_j)
Joining with `by = join_by(county)`
pal <- colorRampPalette(c("#ff1a1a", "#ffa142", "#ccff66"))
mapview(life_exp_md,
        label=life_exp_md$county,
        zcol="Life_Expectancy",
        col.regions=pal(6),
        at = (seq(min(life_exp_md$Life_Expectancy), max(life_exp_md$Life_Expectancy), length.out = 6)),
        legend = TRUE,
        layer.name = "2018-2020 Life Expectancy in Years",
        popup=popupTable(life_exp_md, zcol = c("Life_Expectancy")))